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ABSTRACT 

With the launch of several planetary missions in the last 
decade, a large amount of planetary images has already 
been acquired and much more will be available for anal- 
ysis in the coming years. The image data need to be ana- 
lyzed, preferably by automatic processing techniques be- 
cause of the huge amount of data. Although many au- 
tomatic feature extraction methods have been proposed 
and utilized for Earth remote sensing images, these meth- 
ods are not always applicable to planetary data that often 
present low contrast and uneven illumination character- 
istics. Different methods have already been presented for 
crater extraction from planetary images, but the detec- 
tion of other types of planetary features has not been ad- 
dressed yet. Here, we propose a new unsupervised method 
for the extraction of different features from the surface 
of the analyzed planet, based on the combination of sev- 
eral image processing techniques, including a watershed 
segmentation and the generalized Hough Transform. The 
method has many applications, among which image regis- 
tration and can be applied to arbitrary planetary images. 

Index Terms — Crater Detection, Feature Extraction, 
Watershed Segmentation, Hough Transform. 

1. INTRODUCTION 

The Lunar Reconnaissance Orbiter (LRO) is a NASA mis- 

sion, aimed at creating a comprehensive atlas of the features 
and resources of the Moon to prepare exploration and scien- 
tific missions to the Moon and to study its environment [1]. 
Having officially reached lunar orbit on June 23rd, 2009, LRO 
has now marked one full year on its mission of scouting the 
Moon. In only the first year of the mission, LRO has gathered 
more digital information than any previous planetary mission 
in history. Different types of data are being collected by LRO 
at different times, by different sensors, and from different 
view-points: Multitemporal, multimodal and stereo-images 
will be soon available. Therefore, image registration will be 
an essential task to jointly exploit, integrate, or compare all 


these different data. Feature extraction, i.e., extraction of spa- 
tial features in the images, is the first step in the image regis- 
tration process. Furthermore, the feature extraction is impor- 
tant for further analysis of the data. 

Identification of spatial features on planetary surfaces can 
be manually performed by human experts but this process can 
be very time consuming. Therefore, a reliable automatic ap- 
proach to detect the position, structure, and dimension of each 
feature is highly desirable. This is a difficult task for several 
reasons: Limited data are usually available, the quality of the 
images is generally low (i.e., it depends on illumination, sur- 
face properties and atmospheric state), and the features that 
are present in the images can be barely visible due to atmo- 
spheric erosion and they may be based on different structure 
types of variable sizes. 

Among the typical features in planet-surface imagery, 
craters play a primary role. Detection of craters has been 
widely addressed and different approaches have recently 
been proposed in the literature, based on the analysis of 
planetary topography data [2], satellite images [3] in die vis- 
ible spectrum and the infrared spectrum. Here, we focus on 
optical image-based approaches for crater and rock detec- 
tion. Different approaches have been presented, based on 
template matching [3], [4], [5], texture analysis [6], neural 
networks [7], [8], boosting approaches [9], or a combina- 
tion of these techniques [10], [11]. In particular, in [12], 
the identification of impact craters was achieved through the 
analysis of the probability volume created as a result of a 
template matching procedure. Such methods enable the iden- 
tification of round spatial features with shadows. Kim and 
Muller [6] presented a crater detection method based on tex- 
ture analysis and ellipse fitting. That method was not robust 
when applied to optical images, hence it was performed by 
using fusion techniques exploiting both DEM and optical 
images. In subsequent work [13], in order to automatically 
detect craters on Mars, the authors proposed a combination 
of edge detection, template matching, and supervised neural 
network-based schemes for the recognition of false positives. 
In a different approach, Martins et al. [14] adopted a super- 



vised boosting algorithm, originally developed by Viola and 
Jones [15] in the context of face detection, to identify craters 
on Mars. In [16], Urbach and Stepinski presented a differ- 
ent approach for crater detection in panchromatic planetary 
images. The method in [16] is based on using mathematical 
morphology for the detection of craters and on supervised 
machine learning techniques to distinguish between objects 
and false alarms. Until now, apart from craters, no other types 
of planetary features have been addressed in the literature. 

In order to overcome the typical problems of planetary im- 
ages with a lack of contrast, poor illumination, and a lack of 
good features, we propose here a new unsupervised region- 
based approach for the extraction of different planetary fea- 
tures. The main contribution of this paper is a novel unsuper- 
vised method for the segmentation of planetary images, aimed 
at extracting not only craters but also other planetary features. 
Moreover, the approach allows not only to locate the features, 
but also to reconstruct their shape. 

The proposed feature extraction approach is described in 
Section 2. Section 3 includes the presentation and analysis 
of the obtained results. In the same section, the application 
of the proposed technique to an image registration process is 
presented. Experimental results are shown with images rep- 
resenting the surface of Mars, because LRO images are not 
available yet. Finally, conclusions and ideas for future im- 
provements are presented in Section 4. 

2. THE PROPOSED APPROACH 

Planetary images show the surface of a planet and its struc- 
tures. The aim here is to automatically detect the different 
structures that are present on a considered planetary surface 
by using image analysis techniques. The extracted features 
could then be used for image registration registration pur- 
poses , as will be shown in Section 3. Different types of fea- 
tures are present in the planetary images, and their size, shape 
and position are estimated by applying different methods. 

The main features to be extracted are craters and rocks. 
Craters are objects of approximately elliptical shape with 
shadows, due to their depth and uneven illumination. Rocks 
have small circular or elliptical shape, with almost no shad- 
ows. The extraction of these spatial features is a difficult task, 
because planetary images are blurry, quite noisy, present lack 
of contrast and uneven illumination, and the represented ob- 
jects are not well defined. For these reasons, a region-based 
approach, based on segmentation, has been chosen in order 
to address such problems. Segmentation is the process of 
partitioning an image into multiple regions, for instance, in 
order to distinguish objects from the background. A frequent 
approach to segmentation introduces a set of characteristic 
points that are related to the objects to be detected, auto- 
matically selected and used as “seed points” to segment the 
images. Many segmentation approaches have been explored 
in the literature. Here, the watershed algorithm, presented 


by Beucher in [17], has been chosen, a method which is 
automatic, robust and fast. The basic concept of watershed 
segmentation is visualizing a grey-level image into its topo- 
graphic representation (i.e., the grey level of a pixel represents 
its elevation). A flooding process starts from the minima of 
the image, so that the merging of the flooding coming from 
different sources is prevented. As a result the image is parti- 
tioned into two different sets: The catchment basins (i.e., the 
regions) and the watershed lines (i.e., the region boundaries). 
The flowchart of the proposed technique for spatial feature 
extraction is shown in Figure 1. 



Fig. 1. Flowchart of the proposed approach. 


Before applying feature extraction techniques, the images 
need to be preprocessed. First, the noise is smoothed by ap- 
plying Gaussian filtering and median filtering operations in 
cascade [18]. Then, in order to detect edges, the image gradi- 
ent is computed by using the Canny edge detector [19]. As an 
intermediate result of this operation an intensity gradient, J 3 , 
is generated. Then, by applying a non-maximum suppression 
algorithm followed by an hysteresis thresholding to a bi- 
nary gradient image, h, is obtained but this image shows the 
contours of the objects represented in the original image. 

Rocks generally appear like closed contours in /&, be- 
cause of the absence of shadows. In order to extract these 
features, the watershed segmentation algorithm is applied to 
lb and closed contours are extracted. All the areas included 
within a closed contour correspond to “seed point-areas,” and 
are identified as regions. The result of this first step is a binary 
image R that shows the rock boundaries. 

While rocks appear like closed contours and can be easily 
detected, craters have a more complex structure and, due to 
their depth and uneven illumination, exhibit shadows. Then- 
borders can be approximated with incomplete non-continuous 
elliptical curves. A generalized Hough accumulator [20] is 
used to identify the seed points to detect these structures from 
7fc. For every pair of pixels that are detected as edge points 
in I b and exhibit opposite gradient directions, an accumula- 
tor, corresponding to the median point between them in the 
image plane, is incremented of a unit value. The maxima of 
the accumulator are taken as centers of ellipses. The three 
parameters describing the ellipse centered in each detected 
maximum are then computed and a 3D accumulator is used 
to estimate the two semi-axes and the direction angle of the 
ellipse from all the pairs of points that contributed to the ac- 
cumulator in the considered center. The center of each ellipse 
that has been generated is used as a seed point for segmenta- 
tion. Starting from all the detected seed points, a watershed 






algorithm is applied to I g and the craters are identified. I g 
is used here because it represents not only the edges but also 
the elevation information. As a result, a binary image C that 
shows the crater boundaries is obtained. A binary image, I c , 
which represents the contours of all detected features, is cre- 
ated. The binary image, I c , shows the boundaries of the fea- 
tures, identifies their location and estimates their shape. 

The proposed technique for feature extraction can be used 
to register image pairs representing the same scene. For reg- 
istration, two binary images need to be extracted from both 
images to be registered and their match can be estimated (in 
Section 3 the proposed approach to image registration is pre- 
sented). 


3. EXPERIMENTAL RESULTS 

Since LRO images are not yet available, experiments were 
carried out using Mars data, collected during the 2001 Mars 
Odyssey mission, by the THermai Emission Imaging System 
(THEMIS), an instrument on board the Mars Odyssey space- 
craft. Such an instrument combines a 5-band visual imaging 
system with a 10-band infrared imaging system [21]. Both 
visible (VIS) and infrared (IR) images, with a resolution of 
18 meters and 100 meters per pixel, respectively, were used 
to test the proposed approach. For the experiments 5 VIS and 
7 IR images were selected. 

Reference data were generated by manually analyzing 
each image of the data set and identifying all the craters and 
rocks that are present. Only objects completely included 
within the images were considered (i.e., objects cut by the 
borders of the image were discarded). No limits were im- 
posed on the minimum dimensions of the features to be 
detected. A quantitative assessment of the obtained results 
by the proposed method was performed using these reference 
data. This was accomplished by comparing the obtained re- 
sults with the labeled features in the correspondent reference 
map, RM. The Detection percentage D, the Branching factor 
B, and the Quality percentage Q were computed as follows: 

100 *TP FP 100 * TP 

~ TP + FN ’ ** ~ TP ’ ^ _ TP + FP + FN ’ 

where True Positive (TP) is the number of detected features 
that correspond to labeled objects in RM , False Positive 
(FP) is the number of features detected by the proposed 
approach, which do not correspond to any object in RM , and 
False Negative (FP) is the number of objects in RM that 
have not been detected by the proposed approach. The global 
values of D, B , and Q and the total number of TP , FP , 
and FN obtained by the proposed approach both for VIS and 
IR data are shown in Table 1. The global values of D for 
VIS data and IR data were about 82% and 78%, respectively; 
these high values indicate a good detection rate (because of 
the high number of TP). B was about 0, 03 for VIS and 0, 05 
for IR, which indicate a small amount of false detections with 


Table 1. Average numerical performance of the proposed ap- 
proach as measured by Detection percentage (D), Branching 
factor (B) and Quality percentage (Q) and total number of 
True Positive (TP), False Positive (FP) and False Negative 
(FN) detections on visible and infrared images, respectively. 


Data 

D 

B 

Q 

TP 

FP 

FN 

VIS 

82% 

0,03 

81% 

197 

6 

58 

IR 

78% 

0,05 

75% 

211 

6 

61 


respect to the true detections in both cases, thanks to the small 
number of FP. The detection performance of the proposed 
approach in terms of D, B, and Q compares favorably with 
most of the results previously published for automatic crater 
detection methods [13, 16, 22]. Nevertheless, no full com- 
parison is possible, since all the other methods are based on 
using only crater features. 

Visual results are shown for a partition of a single band 
VIS image (Figure 2-a). The grey level image is first pre- 



(a) Original image 



(c) “Seed point” watershed 
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(b) Binary watershed 



Fig. 2. Experimental results obtained by applying the pro- 
posed method to a VIS image, (a) Original image, (b) water- 
shed segmentation applied to h, (c) watershed segmentation, 
using the maxima of the Hough accumulator as “seed points”, 
and (d) detected features. Each color in the segmentation map 
denotes a different region. 

processed by a smoothing filter in order to reduce the noise. 
The Canny operator is applied to the smoothed image. Sub- 
sequently, in order to extract the rocks, a watershed algorithm 
is applied to the binary image gradient h- Rock segmentation 
results are shown in Figure 2-b. Then, the generalized Hough 
transform is computed and a watershed segmentation is ap- 
plied, starting the flooding process from the ellipse centers 
and leading to the detection of the craters. The segmentation 






result are shown in Figures 2-c. Finally, the extracted fea- 
tures, including both rocks and craters, are combined into a 
binary map and shown in 2 -d, transparently superimposed to 
the original image. By a visual inspection, it is possible to 
appreciate the accuracy of both the detection and the recon- 
struction of the feature shape. 

Visual results are also shown for a partition of the first 
band of an IR image (Figure 3-a). Figure 3-b shows the seg- 



(a) Original image (b) Binary watershed 



(c) “Seed point” watershed (d) Feature map 


Fig. 3. Experimental results obtained by applying the pro- 
posed method to the first band of an IR THEMIS image. 

mentation results when watershed is applied to I 5 . Figure 3-c 
shows the crater segmentation results. The different extracted 
features are combined and shown in 3-d. In this example, 
not all the features are detected. This is because their con- 
tours were not extracted by the Canny operator. An improve- 
ment of the edge detection step may improve the accuracy of 
the method and it is currently under investigation. Anyway, 
we will show below that the detected features are enough to 
achieve an accurate registration. 

To demonstrate the applicability of the proposed method 
to registration, two different non-registered bands of an IR 
image are used. In order to show the results, the same parti- 
tion of Figure 3-a is used; in particular, the 4 th and 5 th bands 
were selected (Figures 4-a and 4-b, respectively). For both 
band images, craters and rocks are detected and their con- 
tours are extracted and represented in binary contour images, 
I Cj 1 and / Cj2 , as shown in Figures 4-c and 4-d, respectively. 
The rotation and translation between the two bands are visi- 
ble by looking at Figure 4-e, in which the two non-registered 
contour images are superimposed in a false-color representa- 
tion. The contours extracted from the 4 th band image, I Ci u 
are represented in green, whereas the 5 ^-band contours, I c , 2 , 
are shown in red. The registration scheme used in this phase 
was based on a global optimization technique aimed at esti- 
mating the optimum parameters of an image transformation 
model. The contour images, which represent the features of 



(a) Original Reference Image 



(c) Contours of (a) 



(b) Original Input Image 



(d) Contours of (b) 



(e) Superimposed Contours 


(f) Registration Results 


Fig. 4. Registration results for a partition of (a) the A th and 
(b) the 5 th bands of an IR image, (c) and (d) feature con- 
tours extracted from (a) and (b), respectively, (e) Contour 
images superimposed and represented in a false-color com- 
position (i.e., the green plane is (c), the red plane is (d), and 
the blue plane is identically zero), (f) Registration results, by 
using a checkerboard representation. 


the two input image bands, were fed as inputs to an opti- 
mization module. The transformation matrix was to be op- 
timized: Its goodness was evaluated by an objective function 
and its optimization was achieved by applying a genetic algo- 
rithm [23]. After the optimum matrix was estimated, it was 
used to transform and interpolate one of the two bands with 
respect to the other one. The co-registered bands are shown 
in Figure 4-f, by using a checkerboard representation: Each 
tile of the board represents the registered input band and the 
reference band, alternately. The registration accuracy can be 
evaluated by looking at the continuity of the features at the 
borders of the tiles. The visual analysis of Figure 4-f suggests 
that the registration performed very well; craters and ridges 
appear continuous at the borders, i.e., the points of overlap. 

4. CONCLUSIONS 

In this letter a novel unsupervised region-based approach has 
been proposed for automatic detection of spatial features that 
characterize planetary surfaces. The proposed approach has 
been applied to the registration of planetary data. 










The features to be extracted are not as well contrasted nor 
defined as for Earth data. However, small rocks, which are 
not affected by uneven illumination, can easily be detected 
by the proposed approach. Crater detection is more difficult 
than the rock detection, because of their depth and spatial ex- 
tent and, consequently, their contours are often blurry and not 
continuous. Nevertheless, we showed here that their identifi- 
cation can be achieved and the proposed approach provided 
quite accurate results. The accuracy of the detection has been 
assessed by comparison to a manually generated reference 
map. The results in terms of several indices based on true and 
false positives compared favorably to previously proposed ap- 
proaches. Moreover, we showed that the extracted features 
can be used to accurately register pairs of image bands ac- 
quired from the same sensor. The accuracy of the registration 
step is confirmed by visual inspection of the results. 

The proposed approach represents the first important step 
for many applications dealing with all the various data that are 
being collected by LRO and other planetary missions, among 
which image registration and image analysis, with the aim of 
selecting safe landing sites, identifying lunar resources, and 
preparing for subsequent explorations of the Moon and Mars 
by both robots and humans. 

In our future work we plan to integrate the shadow infor- 
mation around the features in order to improve the reliability 
of the edge detection and reduce the false alarms in the con- 
tour map 7 C . The approach will be extended to the registration 
of multisensor images. 
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